What is moving in silica at 1 K? A computer study of the low-temperature anomalies 
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Though the existence of two-level systems (TLS) is widely accepted to explain low temperature 
anomalies in many physical observables, knowledge about their properties is very rare. For silica 
which is one of the prototype glass-forming systems we elucidate the properties of the TLS via 
computer simulations by applying a systematic search algorithm. We get specific information in 
the configuration space, i.e. about relevant energy scales, the absolute number of TLS and electric 
dipole moments. Furthermore important insight about the real-space realization of the TLS can be 
obtained. Comparison with experimental observations is included. 

PACS numbers: 61.43.Fs, 63.50.+X, 



Most kinds of disordered solids show anomalous be- 
havior at very low temperatures (Kelvin regime and be- 
low) as compared to their crystalline counterparts. Many 
of the observed features can be explained by the Stan- 
dard Tunneling Model (STM) [J, and its generaliza- 
tion, which is the Soft-Potential Model 

ill 

The basic 

idea of the STM is to postulate the possibility of localized 
transitions between different configurations, i.e. adjacent 
minima of the potential energy landscape. Such a transi- 
tion can be described by a double- well potential (DWP), 
characterized by an asymmetry A, potential height V 
and distance d between both configurations. From a 
quantum-mechanically perspective at low temperatures 
the system is tunneling between both configurations and 
the DWP is characterized by the lowest two eigenstates. 
If their energy difference E is in the Kelvin regime, these 
DWP may contribute to the low-temperature anomalies. 
Then one may speak of Two Level Systems (TLS). The 
TLS can couple to strain and electric fields and therefore 
show up in observables like thermal conductivity, sound 
absorption and dielectric response |5(. Recently, even 
observations about the interaction of different TLS have 
been reported |E 0, 111 • 

So far it has not been possible to derive a theory of 
the low-temperature anomalies of real glass-forming sys- 
tems from first principles, except for mean-field models 
and random first order transition theory |lfj|. Thus 
for a prototype system like SiC>2 the STM has to be ba- 
sically considered as phenomenological. Important ques- 
tions emerge. (Ql) How many TLS are present? This 
is, of course, a central experimental observables. (Q2) 
What are the energetic properties of the TLS, e.g. typi- 
cal barrier heights? (Q3) What is the average microscopic 
nature of the TLS? 

Computer simulations may help to shed some light on 
the nature of the low-temperature anomalies. In pre- 
vious work on Si02 the trajectories, generated either 
by molecular dynamics [III Il2j| or by the activation- 
relaxation technique |13| . have been analyzed with re- 
spect to transition events between different structures. 
Both approaches yield some interesting insight into the 
nature of relaxation processes in SiC"2. (Ql), however, 
requires a systematic search procedure, and (Q2) and 



(Q3) a sufficiently large number of characteristic DWP 
and thus an efficient search method. This has not been 
the scope of previous simulations on Si02. 

In recent years we have developed a set of simulation 
techniques which allow us to approach these questions 
[T3. fl5{ . (Ql) Starting from representative low-energy 
structures we systematically search for adjacent minima 
of the potential energy landscape, i.e. DWP. It is pos- 
sible to formulate an intrinsic completeness criterion to 
check whether an existing adjacent minimum has indeed 
been found |l4j . If this criterion is fulfilled we have access 
to the absolute number of TLS. Saddles are determined 
via a robust saddle-search algorithm 16]. (Q2) The typ- 
ical interatomic energy scale and thus the typical range 
of DWP asymmetries is of the order 1 eV. Thus via a 
direct search it is impossible to find a set of DWP with 
asymmetries in the Kelvin regime. In the Soft-Potential 
Model a DWP is parametrized by a quartic polynomial 
Sn=2 w n% n , thereby reproducing the values for A, V and 
d. If not mentioned otherwise, d corresponds to the mass- 
weighted distance, obtained from a straight connection 
of both minima with the transition state. Using this 
parametrization we can first determine distribution func- 
tions Pi(wi) from our set of numerically found DWP and 
then generate an arbitrary large number of DWP with 
the same statistics. Thus it is possible to estimate the 
DWP distribution function p(d, V, A) over a broad range 
of parameters and in particular to get important infor- 
mation about the properties of nearly symmetric DWP, 
i.e. TLS. (Q3) Due to the effectiveness of the search al- 
gorithm we find a sufficiently large number of DWP to 
perform a reasonable statistical analysis with respect to 
microscopic properties. All technical aspects and in par- 
ticular the justification of the parametrization approach 
can be found in 0, [^. So far we have successfully 
applied these methods to Lennard-Jones model systems 
El 13 El 13 El; see also [H and chapter 8 of Ref.0. 

In this paper these techniques are for the first time 
applied to pure SiC"2, modelled by the BKS-potential ppf . 
Despite the enormous numerical effort, involved in these 
calculations, we were able to obtain detailed results about 
the central questions (Q1)-(Q3), introduced above. 

We have analyzed system sizes of 150 and 600 particles 
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Figure 1: The two structures correspond to the minima of 
a DWP. The black particle is the one, moving most, and the 
white particles mark the 11 most displaced ones, of which two 
oxygens are not shown as they are not connected to the above 
network. The grey particles correspond to the less moving 
oxygen (small spheres) and silicon (large spheres) atoms. 



with the standard density of 2.3 g/cm 3 and a short-range 
cutoff for the BKS-potential of 8.5 A. The starting config- 
urations for the systematic search correspond to equilib- 
rium configurations at 3000 K, which subsequently were 
minimized. We have obtained starting structures with 
and without defects, i.e. deviations from a perfect tetra- 
hedral coordination. It turned out that the properties 
of DWP are totally different in both cases because a de- 
fect is very likely related to DWP with high barriers 0] 
and a stronger spatial localization. Experimental silica 
samples, however, freezing in at the calorimetric T g , do 
not show a high number of defects [2l|. Therefore we 
have only considered structures without defects. Actu- 
ally, these non-defect structures are very similar to those 
one would obtain from simulations at much lower tem- 
peratures 22]. Furthermore we have only taken into 
account DWP with |A| < 1500 K. This already cor- 
responds to the symmetric side of the asymmetry dis- 
tribution 0]. With this constraint we have been able 
to produce a set of 250 DWP for 150 particles and 50 
DWP for 600 particles. One example is shown in Fig. 
^ On average we have found one DWP per 14 defect- 
free starting structures. The completeness criterion is 
fulfilled. To a good approximation the distribution of 
asymmetries is constant (data not shown). Thus the 
number of DWPs with |A| < 1 K can be estimated via 
1/(14 • fc B 1500 K • L 3 ) « 5 • 10 44 /(Jm 3 ) where L is the 
length of the simulation box. 

For the 600 particle system the search procedure is 
complicated by the occurrence of independent DWP, 
which makes a systematic search much more time- 
consuming. Due to the better statistics we report re- 
sults for the 150 particle system. We checked, however, 
that apart from minor variations in the participation ra- 
tio (see below) all properties, discussed in this work, are 
identical for both system sizes within statistical uncer- 
tainty. This is compatible with our previous result that 
the thermodynamics of BKS-Si02 with only 99 particles 
is, apart from trivial sca ling , basically identical to that 
of a macroscopic system |22). 

In Fig.|2K we show the distribution of barrier heights. 
Using the parametrization method we obtain the full dis- 
tribution of DWP, i.e. p(d,V,A). Interestingly, the typ- 
ical barrier heights agrees well with the value of approx. 
500 K, estimated from Brillouin scattering experiments 




Figure 2: a) Distribution of barrier heights for the simu- 
lated data with asymmetries up to 1500 K. Furthermore the 
distribution, generated from our statistical analysis, is shown 
(see text) . b) The calculated distribution of barrier heights for 
TLS. In both cases we also show the average distance between 
two minima for a given barrier height. 



|23j . Using the same constraint | A| < 1500 K the original 
distribution is reproduced; see also Fig.|2K- Next we have 
estimated the tunneling matrix element for every DWP 
in the full distribution p(d, V, A), using the WKB approx- 
imation [24| , and calculated the energy splitting E of the 
two lowest eigenfunctions and the relaxation rate r [24| . 
The subset of DWP with E < 2 K and r < 5s, which 
from now on we denote TLS, is relevant in the Kelvin 
regime. The distribution of barrier heights of the TLS is 
presented in Fig-Et*- This distribution shows a maximum 
around 50 K. The average distance between both minima 
varies with increasing potential height. More generally, 
one observes that very nearby minima have smaller asym- 
metries and barrier heights. This correlation has been 
also observed for the Lennard- Jones system |14l [l5j . 

The effective density of states P e ff is accessible, i.e, by 
sound absorption experiments |25|. Based on our distri- 
bution p(d, V, A) and our systematic search procedure we 
can estimate this value from first principles. Again, for 
the necessary technicalities we have to refer the reader 
to the corresponding analysis of Lennard-Jones systems 
[pil [pel . We obtain for pure defect free silica that 



P oS (E) = 1.3 ■ 10 44 (£/1 K) u u 7(Jm 



(1) 



which is of the same order as the value from the sim- 
ple estimation, reported above. P e a(E « 1 K) can be 
directly compared to experimental |26j values between 
(4 — 7) ■ 10 44 /(Jm 3 ) for Suprasil W (205 ppm impurities) 
and 8 • 10 44 /(Jm 3 ) for Suprasil 1(1250 ppm impurities). 
The dependence on energy is weak but significant and 
can be regarded as a correction to the STM. Although 
for a first principle analysis the agreement of P e g(E = 1 
K) with the experimental data is already remarkable the 
remaining differences can be (at least partly) related to 
different aspects: (i) Using a Frank-Condon factor for 
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the multidimensional renormalization and taking into 
account the entropic effects, obtained from a harmonic 
analysis of the transition state and the minima, may in- 
crease the presented value by a factor of 1.5 according to 
our estimation, (ii) The remaining defects, present even 
at T g , may have some impact on this number. Since for 
pure SiC>2 the number of defects is only roughly known it 
is difficult to estimate this contribution |2l| . Since, how- 
ever, defects are very efficient in forming DWP (data not 
shown) even a concentration of 100 ppm contributes to 
P c g. (iii) The impurities in real silica systems contribute 
additional extrinsic DWP. 

In the remaining part of this paper we would like to 
discuss the real-space properties of the DWP found in 
our simulation. For the DWP in Fig. ^ the tetrahedra, 
involved in the main displacement, basically form a ID- 
chain. Visual inspection of about 50 DWP indicates that 
this reduced dimensionality of the transition is a generic 
feature. One can also see that the number of partici- 
pating particles is rather small. Interestingly, only two 
oxygens per tetrahedron move significantly. In what fol- 
lows we will present a statistical analysis of all DWP to 
extract the average behavior. 

For the average participation ratio different definitions 
can be used 0. We obtain (both for N=150 and N=600) 
(d 2 /dl ax ) » l/(dl ax /d^ a 9, <d7E,d?> * 2 ^ N = 
150) and (dVEjdf) ~ 31 ( N = 600 )- For this specific 
analysis d describes the euclidian distance between two 
configurations, d, is the displacement of particle i and 
dmax the displacement of the most displaced particle. Ac- 
tually, very similar values are reported in 0] an d in H3 
where the tetrahedra relaxation upon applied pressure 
has been monitored. The most displaced particle is oxy- 
gen in all analyzed cases (which is further denoted as cen- 
tral oxygen) and the 4 most displaced particles are exclu- 
sively oxygen in over 99% of all analyzed cases. Averaged 
over all particles, the mean square displacement of oxy- 
gen is 2.3 times larger than that of silicon. Furthermore, 
we have checked that the participation ratio decreases 
with decrea sing distance in analogy to the Lennard-Jones 
system 0, One can estimate that the above val- 

ues for the participation ratio decrease by approx. 25% 
when evaluated for distances of d » 1.2 A for typical 
TLS rather than d « 2.8 A for typical DWP (see Fig. EJl . 
Furthermore, the TLS have an electric dipole moment 
which is also accessible experimentally. Using the partial 
charges of the BKS potential we obtain an average dipole 
moment of 0.65 D which is in excellent agreement with 
the experimental value of 0.5 D |28| . 

In the next step we have analyzed the average local en- 
vironment of the central oxygen. It turns out that before 
and after the flip its average Si-O-Si angle is (148.5 ± 1)° 
whereas at the transition state it is approx. 160°. Thus 
to a good approximation the Si-O-Si bond performs a flip 
rather than a simple rotation. Actually, the average Si- 
O-Si angle for all bonds is 152° in the BKS-system [29j . 
which indicates a structural anomaly around the central 
oxygen. The observed strong dependence of the density 
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Figure 3: Special radial distribution function for the central 
particles of the DWP in comparison to a general g(r). No 
differences can be seen beyond the first peak. 



of states on the inter-tetrahedral structure [30J fits in this 
picture, since the relevant Si-O-Si angle is the most im- 
portant parameter for the inter-tetrahedral structure. 

Furthermore, also the O-Si-0 angles vary during the 
transition. The most displaced tetrahedra shows varia- 
tions of the O-Si-0 angles of 2.7° which is nearly half 
of the overall equilibrium variation of intra-tetrahedral 
angles; see also 27]. However, for this quantity no struc- 
tural anomalies could be observed. One may speculate 
that these distortions are an important process of damp- 
ening the motion and thus to localize the DWP. Another 
structural anomaly can be seen from comparing the ra- 
dial distribution functions around the central oxygen at 
the transition and two minima positions with that of an 
average oxygen; see Fig- El First, we observe a shortening 
of the Si-0 bond during the transition. Thus the Si-O- 
Si bond flip is achieved by a squeezing of the oxygen 
through the two adjacent Si atoms. Second, and even 
more important, the Si-0 bond lengths for the typical 
DWP minima is significantly smaller than observed for 
the bulk, and their distribution is narrower. 

In order to minimize the energy variation during the 
transition it is essential that the total amount of par- 
ticle displacement is as small as possible. Thus tetra- 
hedra rotations around the threefold axis, giving three 
equally large oxygen displacements, should be very un- 
favorable. A better choice would be rotation around one 
of the tetrahedral edges, thus moving only two oxygen 
atoms. This type of motion has already been indicated 
in Fig. ^ where only two oxygens at most tetrahedra are 
largely displaced. 

To establish a statistical picture of the typical tetrahe- 
dron rotation the average displacements of the 4 oxygen 
atoms at each tetrahedron have been investigated. In 
Fig. 01 results are presented for the 10 most displaced 
tetrahedra averaged over all TLS. To give a better im- 
pression of the rotation the three differences of the dis- 
placements are shown. The first dot in each section corre- 
sponds to the difference of the displacements of the two 
most displaced oxygens and so on. In agreement with 
expectation the nature of the rotations is quite different 
to case A (rotation around threefold axis). Interestingly, 
the most displaced tetrahedron indeed shows a rotation 
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Figure 4: DifFerences of the displacements of the four oxygens 
at the 10 most displaced tetrahedra averaged over all DWP. 
Within each tetrahedron the left point displays the displace- 
ment difference (Ado,i) between the most displaced oxygen 
and the second most displaced oxygen. The middle point 
represents Adi t 2 and the right point Acta, 3- The leftmost (A) 
sequence of points corresponds to a threefold rotation. The B 
sequence corresponds to a rotation around a tetrahedra edge. 

with some similarity to the proposed rotation axis along 
a tetrahedral edge (case B in Fig. 0J. All other rota- 
tions follow a different scheme, which involves four dif- 
ferent displacements for the oxygen atoms. This can be 
possibly interpreted as a transition to a statistical (low- 
amplitude) displacement with no symmetry criterion for 
the rotation axis. 

One can compare the observed rotation scheme with 
the proposed motion of five connected tetrahedra by 
Buchenau et al. which is often used as a microscopic 
visualization for Si02-TLS [3l|]. The different displace- 



ments for all oxygens connected to a tetrahedron as well 
as the Si-O-Si flip motion fully agree with our results. 
There are, however, some important modifications sug- 
gested by this work: (i) the individual tetrahedron ro- 
tation axes does not always go through an oxygen or 
silicon atom, (ii) tetrahedron distortion should be con- 
sidered, (iii) the spatial structures are one dimensional 
rather than three dimensional. It might be also interest- 
ing to compare these modes in detail with the localized 
vibrational modes studied, e.g., in [3^ . 

In summary, for the first time it is possible to ob- 
tain a systematic numerical description of the TLS for 
BKS-Si0 2 . The density of TLS and the electric dipole 
moment are in semi-quantitative agreement with exper- 
imental data. Important microscopic properties of the 
TLS like structural anomalies and the nature of the dis- 
placements have been gained. In principle, microscopic 
information about TLS is also available from recently 
performed experiments on magnetic field dependent po- 
larization echoes. Thus, with the complementary views 
from experiment and simulation a detailed microscopic 
picture of the low-temperature anomalies will hopefully 
emerge. 
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